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Abstract It has long been suspected that the Global Navigation Satellite Systems exist in a background of 
complex resonances and chaotic motion; yet, the precise dynamical character of these phenomena remains 
elusive. Recent studies have shown that the occurrence and nature of the resonances driving these dynamics 
depend chiefly on the frequencies of nodal and apsidal precession and the rate of regression of the Moon’s 
nodes. Woven throughout the inclination and eccentricity phase space is an exceedingly complicated web¬ 
like structure of lunisolar secular resonances, which become particularly dense near the inclinations of the 
navigation satellite orbits. A clear picture of the physical significance of these resonances is of considerable 
practical interest for the design of disposal strategies for the four constellations. Here we present analytical 
and semi-analytical models that accurately reflect the true nature of the resonant interactions, and trace 
the topological organization of the manifolds on which the chaotic motions take place. We present an atlas 
of FLI stability maps, showing the extent of the chaotic regions of the phase space, computed through a 
hierarchy of more realistic, and more complicated, models, and compare the chaotic zones in these charts 
with the analytical estimation of the width of the chaotic layers from the heuristic Chirikov resonance- 
overlap criterion. As the semi-major axis of the satellite is receding, we observe a transition from stable 
Nekhoroshev-like structures at three Earth radii, where regular orbits dominate, to a Chirikov regime where 
resonances overlap at five Earth radii. From a numerical estimation of the Lyapunov times, we find that 
many of the inclined, nearly circular orbits of the navigation satellites are strongly chaotic and that their 
dynamics are unpredictable on decadal timescales. 

Keywords Medium-Earth orbits • Secular dynamics • Orbital resonances • Chaos • Fast Lyapunov 
Indicators • Stability maps 


1 Introduction 


Among the earliest problems in astrodynamics were, quite naturally, those concerned with the motion of 
an artificial satellite in the gravitational field of the oblate Earth; and indeed their solutions in the hands 
of Brouwer, Garfinkel, Kozai, and others were intimately bound with the rigorous development of artificial 
satellite theory (Jup pjl988 l. These considerations placed emphasis on the construction of increasingly more 
accurate analytical theories, in the style of our forebears, valid in idealized situations or on short mission 
timescales. Such an approach was of course justified by the need of astrodynamical practice, but in recent 
years astrodynamics has had to face new problems concerning the long-term motion of space debris which 
forced it to abandon the somewhat utilitarian investigations mentioned above. This renaissance has been 
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spurred on not only by the sobering implications of space debris, but also by the ability to perform computer 
simulations of ever-greater sophistication and by the realization that chaos has played a fundamental role 
in the dynamical evolution of the Solar System (Morbidelli 2002). 

The proliferation of space debris orbiting Earth has stimulated a deeper understanding of the dynamical 
environments occupied by artificial satellites (B reiter| 2001), and the smaller disturbing forces arising from 
lunisolar gravity have became particularly interesting now (Breiter 1999). For near-Earth satellites, the 
perturbing effects of the Sun and Moon, often assumed negligible in comparison to that of the Earth’s 
oblateness, may have small consequences early but profound consequences late. The subtle effects that 
accumulate over long periods of time are known as resonances and while their basic theory was developed 
during the early Sputnik era (Upton et al. 1959|), a detailed and systematic investigation of their dynamical 
effects is hitherto missing. The problem is particularly timely since, in the past decade, a number of 
numerical studies have shown that the medium-Earth orbits (MEOs) and graveyard regions of the navigation 
satellites are unstable (q.v. D eleflie et al.|poTT and references therein). This instability manifests itself as 
an apparent chaotic growth of the eccentricity on decadal to centennial timescales. 

To identify the underlying dynamical mechanism responsible for the noted instability in the MEO region, 


Rosengren et al. (2015) investigated the main resonances that organize and control the long-term dynamics 


of navigation satellites. They confirmed the complex resonant structure identified nearly twenty years ago 
by Ely and Howell (1997), which has been unjustifiably overlooked by the space debris community, and 
showed that these lunisolar resonant harmonics together with those studied by Delhaise and Morbidelli 
(1993) interact to produce a dense stochastic web in inclination and eccentricity phase space. The close 
spacing of the resonant harmonics suggested the potential for significant chaos in the orbital eccentricities 
and inclinations. This was partially demonstrated numerically through stroboscopic Poincare maps, showing 
these resonances to be the preferential routes for transport in phase space, by which orbits can explore large 
regions jumping from one resonance domain to another ( Rosengren et ah| 2015l. The nature of the chaotic 
layer, however, was not explored. 

The present investigation is a first attempt towards understanding the chaotic structure of the phase 
space near lunisolar secular resonances. This is studied analytically by examining the width of the chaotic 
layer from the overlap of neighboring resonances ( Chirikov|1979 Mardling 2008) and numerically using the 
fast Lyapunov indicator (FLI) (Froeschle et al.|1997). To cast light on the speculations made in Rosengren 


et al. (2015) and to gain insight into the structure, extent, and evolution of the chaotic regions, we compute 


FLI stability maps using a hierarchy of more realistic dynamical models. We make a detailed comparison 
between the analytical and numerical results, and test the predictions of the heuristic Chirikov criterion in 
systems of more than two degrees of freedom. 


2 The Chirikov resonance-overlap criterion 


The long-term evolution of resonant orbits can be quite complex, and what has become clear over the past 


few decades is that the phenomenon of chaos is tied fundamentally to the dynamics of resonances (Chirikov 


1979 Morbidelli 2002). As an instance we can mention the problem of motion of asteroids in the main 


belt between Mars and Jupiter on which so much of our astronomical knowledge is based; and here the 
light thrown by purely dynamical considerations upon, for example, the problem of Kirkwood’s gaps in the 
distribution of the asteroidal semi-major axes reveals the sorts of resonant interactions that have helped to 
shape the Solar System. The physical basis of chaos and instability in nearly integrable, multidimensional 
Hamiltonian systems, of which the Sun-Jupiter-Saturn-Asteroid four-body problem is one example, is the 
overlapping of nonlinear resonances (Chirikov 1979). 

The overlapping of resonances as the mechanism of onset and evolution of chaos has been known since 
the late 1960s. This important result from nonlinear dynamics can be made intelligible without entering 
into the lengthy and subtle mathematical discussion demanded by an exhaustive treatment of the subject 
(Mardling][2008). Astronomers long ago noted the connection between the dynamics of a single (isolated) 
resonance and the libration of a mathematical pendulum. An expansion of the Hamiltonian around the 
resonant location and a canonical transformation reduces the Hamiltonian to a pendulum-like system. 
The resonance angle librates when the frequencies are nearly commensurate; the domain or width of the 
resonance being the distance from exact resonance to the separatrix (a homoclinic curve separating the 
domain of circulation and libration). Whereas a single resonance exhibits a well-behaved, pendulum-like 
motion, when the system is moving within a region of the phase space where two or more resonances are 
present, the single resonance theory breaks down and any attempt to describe the interaction rigorously 
is bound to fail. The nature of this breakdown manifests itself in numerical solutions as zones of chaotic 
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motion near the boundary of each interacting resonance. The Chirikov resonance-overlap criterion simply 
states that chaos will ensue if two dynamically relevant harmonic angles, when neglecting the interaction 


between them, are both analytically calculated to be librating in the same region of phase space (Chirikov 


1979). 


2.1 Location of lunisolar resonance centers 


Our aim here, following Ely and Howell (19971 and Rosengren et al. (2015), is to identify the significant 
resonances so that we can apply the resonance-overlap criterion and analytically determine stability bound¬ 
aries. We focus on the MEO region between roughly three and five Earth radii and, as a result, are permitted 
to make certain approximations and assumptions. There are two principal classes of resonances which affect 
the motion of satellites in these orbital altitudes: tesseral resonances, where the orbital periods (or mean 


motions) of the satellites are commensurable to the Earth’s rotation rate (Lane 


1988 


Celletti and Gale§ 


2014), and lunisolar secular resonances, which involve commensurabilities amongst the slow frequencies of 


orbital precession of a satellite and the perturbing body (Ely and Howell 1997). The dynamical system 


with tesseral, zonal, and lunisolar perturbations is, however, characterized by two different timescales and 
two independent modes. Consequently, we focus here only on the resonant effects of secular origin, which 
dominate the eccentricity and inclination evolution on long timescales (Batygin et al. [20151. 

The analytical methodology for computing the lunisolar resonant effects is based on a theoretical series 
development of the lunar and solar disturbing functions (the negative potential functions of the perturbing 
accelerations). These harmonic series expansions are expressed most simply, in their dependence on the 
orbital elements, when the latter are defined with respect to the ecliptic plane for the Moon and with 


respect to the equatorial plane for the satellite and the Sun (Hughes 1980). For our purposes, we can 


truncate these series to second order in the ratio of semi-major axes, so that the lunar and solar potentials 
are approximated by quadrupole fields. In the secular approximation, short-periodic terms which depend 
on the fast orbital phases (i.e., the mean anomalies of both the satellite and the perturbing bodies) can be 
readily averaged over and dropped from the disturbing functions ( Arnold et aL| 2006). 

The secular and quadrupole order of approximation for the lunar disturbing function expansion follows 


from Lane (1989) and has the form 


with harmonic angle 


TZm = EEE ^2—2p,m,±s COS 2p,m,±s? 

m =0 s=0 p =0 


@2-2p,m,±s = (2 - 2p)w + mQ ± s(Qm - tt/2) - y s ir, 


( 1 ) 

( 2 ) 


and associated harmonic coefficient 
2 


h 


M 


2—2p,ra,±s — 


fl M CL 


r(—1) 


|m/2j (2 ~ S)! 


<(l-4) 3/2V ^ £m (2 + m)! J 

where m, s, and p are integers, [ J is the floor operator, and the quantities t m and y s are such that 

1 if m = 0, 


F2,mA i ) F 2,sAiM)H 2 ^ 2p -2(e)(-l) m+s U™’ Ts (e), (3) 


Cm — 


2 if m ^ 0, 


f 0 if s is even, 

1/2 if s is odd. 

The semi-major axis a, eccentricity e, inclination i, argument of perigee lo, and longitude of ascending 
node 12 are the satellite’s orbital elements relative to the Earth’s equator. The subscript M refers to the 
Moon’s orbital parameters (*m and 12 m are measured relative to the ecliptic plane), /tm is the Newtonian 
gravitational constant times the Moon’s mass, and e is the obliquity of the ecliptic. The quantities T 2 ,m,p(*) 
and F 2 ,s,i (*m) are the Kaula inclination functions (Kaufa 1966), H 2 .p, 2 p- 2 (e) is the zero-order Hansen 
coefficient Xg’ 2 ~ 2p (e) (Hughes 19801, and U™’ s is the Giacaglia function, required for the mixed-reference- 
frame formalism (Giacaglia 1974): 


/_i\m=Fs _i2is r 

t 7-m,ihs v 1 J / _ _ _ /o\m±s / ■ /<-> \±s—ra u I ry2—m / r? -i \ 24-m I 

U ‘ 1 = (2±a)! ( C0S£ / 2 ) (sme/2) p^ 2 ±s \ Z ( Z “ X ) ) 


( 4 ) 
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where Z = cos 2 e/2. Note that the expansion ([I]) is valid for arbitrary inclination and for any eccentricity 
e < 1, and, moreover, this formulation allows us to consider the lunar inclination as a constant and the 
longitude of the ascending node of the Moon’s orbit as a linear function of time (qq.v. Lane 1989 iHughes 


1980). 


The form of the solar disturbing function is considerably simpler, with the secular and quadrupole order 
of approximation following from Kaula’s classical spherical harmonic expansion: 


2 2 

K S = EE ^2—2p,m COS ^2— 2p,mi 

m =0 p =0 


with harmonic angle 

^ 2—2 p,m = (2 - 2p)u 1 + m{n - n s ), 

and associated harmonic coefficient 


h 


s 

2—2 p,m 




z|(i — e§) 3/2 (2 + m)! 


(2 — mV. 

CmTTr:- nF 2 ,m,p(i)F 2 ,m,l(is)H 2 ,p, 2 p- 2 (e). 


(5) 

( 6 ) 

(7) 


Here equatorial elements are used for both the satellite and the Sun. The apparent orbit of the Sun can be 
considered, for all practical purposes, as Keplerian. 

A lunar secular resonance occurs when a specific linear combination of the averaged precession frequen¬ 
cies of the angles appearing Eq. [T]is zero. That is, 

V'2-2 P ,m,±s = (2 - 2p)cj-b m/2 ± sl2 M « 0. (8) 


For the Moon, the rate of change of 17 m, due to the perturbations of the Sun, is approximately —0.053 
deg/d. Since b?s = 0, it follows that secular resonances of the solar origin are characterized by the simpler 
relation 


(2 — 2 p)lu + mfi « 0 . 


(9) 


Under these approximations, the occurrence and nature of the secular resonances driving the dynamics 
depend only on the frequencies of nodal and apsidal precession and the rate of regression of the Moon’s 
nodes. The Moon’s argument of perigee does not explicitly appear in the expansion: this important result 
does not require the additional assumption made by previous researchers ( Delhaise and Morbidelli] 1993| 
Ely and Howell 19971 that the Moon’s orbit be circular. 

For the medium-Earth orbits considered here (i.e., semi-major axes between three and five Earth radii), 
the oblateness effect is at least an order of magnitude larger than that of lunisolar gravity (Ely and Howell) 
1997). We are therefore permitted in the first approximation to neglect the changes in w and 17 from 


lunisolar perturbations in order to make a more efficient mathematical treatment of the problem possible 
without at the same time deviating too drastically from the actual conditions prevailing in the system. The 
secular rates of change of a satellite’s argument of perigee and longitude of ascending node cased by the J 2 
harmonic are such that (Cook|[l9621 


3 T / R\ 5 cos 2 i — 1 

w = I j2n UJ (1-eT ’ 

3 f R\ 2 cos i 

J ~~2 J2n {a) (1 — e 2 ) 2 ’ 


( 10 ) 


where R is the mean equatorial radius of the Earth and n is the satellite’s mean motion. As the semi-major 
axis is constant after averaging, Eqs. m define analytical curves of lunisolar secular resonances in the 
inclination and eccentricity phase space (Fig. [l|. It is particularly noteworthy that the solar resonances and 
the lunar resonance with s = 0, both occurring at 46.4°, 56.1°, 63.4°, 69.0°, 73.2° and 90°, are independent of 
the orbit eccentricity ( Cook|1962 Hughes 1980J; these inclination-dependent-only resonances are intimately 
related to the critical inclination problem (q.v. Jupp 1988). Each of the critical inclinations split into 
a multiplet-like structure of lunar resonance curves corresponding to s € {1,2}, emanating from unity 
eccentricity. As the semi-major axis of the satellite is receding from three to five Earth radii, the resonance 
curves began to intersect indicating locations of multiple resonances, where two or more critical arguments 
have vanishing frequencies. The complex web of intersecting lines suggests the potential for chaos and large 
excursions in the eccentricity, on account of the resonance-overlap criterion. 
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2.2 A measure of the resonant islands 


Each of the resonances in Fig. [l| determines its own domain in the phase space—the resonance width 


Garfinkel 

1966 

Lane 

1988 

Breiter 

CO 

0 

0 

CN 


Breiter 2003). Based on the isolated resonance hypothesis, we derive here a 


measure of the libration islands, expressed in terms of the Delaunay and Keplerian variables. The former 
set of canonical variables are the most frequently employed in astronomical perturbation methods, and 
allow us to appeal to the Hamiltonian formalism. The conjugate coordinates are the Keplerian angles M , 
to, and 17 and the conjugate momenta are functions of the orbital elements a, e, and i, as follows: 


L = y/Jia, l = M , 

G = LV 1-e 2 , g = co, (11) 

H = Gcosi, h = i 7, 

in which p denotes the Earth’s gravitational parameter. 

We take the averaged Hamiltonian H to include the second zonal harmonic of the Earth’s gravitational 
potential and the lowest-order term, i.e., the second harmonic, in the expansion of each of the solar and 
lunar disturbing functions: 


H(G, H , g, h, i; L ) — Wkep + 'H-3 2 + Wm + "Hs, 


( 12 ) 





Fig. 1: The location of resonance centers of the form ip 2 - 2 p,m,±s = (2 — 2 p)io + mfi ± sQm = 0, where 


only the effects of the J 2 perturbation on 10 and 17 have been considered (adapted from Rosengren et al. 


2015). These resonances form the dynamical backbone of the phase space, organizing and controlling the 


long-term orbital motion of MEO satellites. 
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where "Hkep = — gr /2L? is the Kepler Hamiltonian, which depends only on the mean semi-major axis, treated 
as a parameter (M is a cyclic variable and consequently a is constant), Hj 2 is the classical Hamiltonian for 
the averaged oblateness perturbation (Breiter|j2001) 




J 2 R V G 2 - 3 H 2 
~ 4 1A G 3 ’ 


(13) 


Hm(G, H, g, h, t; L ) = —72-m, and Hs(G, H, g, h, t: L ) = — TZs, given by Eqs.[l]and[5j respectively. Equation [l2| 
is a non-autonomous Hamiltonian of two degrees of freedom, depending periodically on time through the 
Moon’s nodal motion. In order to render the Hamiltonian autonomous, we extend the phase space to three 
degrees of freedom by introducing the canonical variables (T, r), such that f = dT-L/dF = %, where the 
new Hamiltonian (without the Keplerian part) takes the form 


H(G,H,r,g,h,r;L) 


H sec (G, H- L ) + -H lp (G, H, g, h, r; L) + fi M F, 


(14) 


in which 


H sec (G, H- L ) = Hj 2 + n s ^ c + HT, (15) 

H lp (G,H,g,h,r,L)=H l & + H 1 i. (16) 

Secular perturbations due to the Moon and the Sun are related to terms with 2 — 2p = 0, m = 0, and s = 0; 
that is, ip 2 - 2 p,m.±s = 0. Long-periodic lunisolar perturbations are connected with ib 2 — r 2 y.m,±s i=- 0. 

Introducing the frequency vector 0 = ( g,h,r ), the lunisolar resonances, Eqs. |8jand|9j may be summa¬ 
rized by the existence of a vector n £ Z* such that 

n ■ 0 = mg + n 2 h + ri 3 T ss 0. (17) 

In order to study each resonance, we treat them in isolation and introduce new action-angle variables, 
appropriate to the resonance involved, by means of the canonical unimodular transformation (linear trans¬ 
formation belonging to SL(3,Z)). The non-null vector n has at most two zero components, and depending 
on the case where m = 0 or n 2 = 0 , we can introduce the resonance angle a = mg + n 2 h + through the 
following unimodular transformations: 

( m n 2 ?i 3 \ / 1 0 0 \ 

0 nj -1 0 , T 2 = I m n 2 713 ■ (18) 

001 / yOOn" 1 / 

The action of T., on 0 gives (T, ■ @ T ) ; = a, i € {1, 2}, where (‘X, ■ 0 T ) ; denotes the i-th components of the 
vector Xj ■ 0 T £ R 3 . To keep the system canonical, new actions A £ R 3 are introduced as 

H = X- T .(G,H,r) T , (19) 


where X, T denotes the inverse transpose of X. ( . These two transformations lead to the new set ( A,<t ) of 
action-angle variables. Using Xi, we have 


A i = n)” 1 G, 

A 2 = —n 2 G + n\H, 

A3 = —n^ 1 ri3G + r, 


ai = a = mg + n 2 h + U 3 T, 
m = n'^ 1 h, 

0-3 = r, 


( 20 ) 


and using X 2 gives 


Ai = G-nm" 1 #, 
A 2 = n 2 1 F[, 

A 3 = —ri 3 FI + n 2 F, 


<n = 9, 

u- 2 = a = nig + n 2 h T n$T, 
13 = 


( 21 ) 


In either case, under the isolated resonance hypothesis, the Hamiltonian, Eq. |14[ expressed in terms of the 
new variables (20 or 21 ) and neglecting the constant terms, reduces to a single degree of freedom with the 


resonant angle a and its conjugated action. Denoting these by (X,x = a), the new Hamiltonian has the 
simple form 


H = f 0 (x)+f 1 (x) cos (x + p). 


(22) 
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with 


r fo(X) = H sec (X) 
\/i(X) =-hn{X), 


where p is a constant phase, H sec is given by Eq. |15[ and hn is the harmonic coefficient in the lunar and 
solar disturbing function expansions, associated with the harmonic angle which is in resonance. For the 
secular resonances with nz ^ 0, hn = h^_ 2v m ±s, as given by Eq. 3 For the lunisolar resonances where 

"3 "0 

(hn) = (h 2 l -2p,m,o) + (h 2 - 2 p,m) + 2(h|L 2 p,m,o) (h|- 2 p,m) cos (mf2s) • (24) 


1966 


Equation |22| is the familiar form for the Hamiltonian when only one critical term is present (Garfinkel 


Lane||l988 BreiterpOOfS). An expansion of the Hamiltonian around the resonant location X* further 


reduces it to a pendulum-like system, which provides the calculation of the amplitude of the libration region 
around the resonance]^] Expanding H in a Taylor series about X*, neglecting terms of 0((X — X*) 3 ) and 
dividing through by —d 2 'H aec /dX 2 1 , gives 


H{I,x) 



+ V cos (x + p), 


(25) 


where I = X — X* and 


hn(X*) 

d 2 H sec /dX 2 \ x=Xt 


(26) 


It should be noted that this last division changes the physical timescale. In this form, the maximum 
excursion in the action, AI , measured from the line / = 0 to the extreme point of the separatrix (the 
resonance half-width), is given by 


AT = 2-Jv. 


(27) 


The width can be expressed in terms of the Delaunay variables using the appropriate transformation, Ti 
or T 2 . It turns out that both transformations, while formally needed, actually lead to the same formulas, 
and thus, in the following, we detail only one of them. Note that the Hamiltonian H is CT 2 and <73 cyclic, so 
that A 2 and A 3 are constants of motion. From Eq. [27] and Eq |20| we find 


AG = |2ni\/i/| . 

Since A 2 = A 2 ,*, we have H = n^r^G + Taking the first variation gives 

AH = \n^n 2 AG\ = \2n 2 yfr\ . 


(28) 


(29) 


The Keplerian elements are related to the Delaunay variables through Eq. m where we note in particular 


that 

From the first variations, we have 
Ae = 

Ai 


■ 


cos l = 


H 

G' 


2 niy/l - el^/v 


L*e.-k 

? 

2 (ri 2 — ni cos i*) yjv 

L*sini*y / 1 - 

el 


e*zie 


1 -el 


n 2 


n 1 sinj* 


— cot i* 


(30) 

(31) 


Ely and Howell 


(1997 J and is more suitable for the calculation 


1 Note that this refinement was not taken into account in 
of the widths for the inclination-dependent-only resonances, where solar perturbations play a non-negligible role (see 
Appendix |b|). 

2 Recall that we need the resonance width in order to determine when neighboring resonances overlap and hence when 

a system is unstable, a la Chirikov. 
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in agreement with the generating-function approach used by Ely and Howell (1997), for which we recall 
here that m = 2 — 2 p and n 2 = m. 

A few words are in order regarding the calculation of v appearing in the above expressions. To be 
consistent with our computation of the resonant centers in Section [2.1[ we neglect the secular contributions 
from the lunisolar perturbations, so that 


H sec (Ai, A 2 ) — 'Hj 2 (Ai, A 2 ) 


J2AV 

4 L 3 n? 


|^(1 — 3% 2 n|)/l 1 3 — 6 n 1 3 n 2 A 2 A 1 4 — 3n 1 4 A 2 A 1 . 


Since X = A \, we finc0 


d 2 U aec 

dX 2 


3 J 2 f?V 

2 L 3 G 5 

3 J 2 R z 


n\ (2-15^3 j +10mn 2 ^-n| 


2a 4 (1 - e 2 ) 5 / 2 


j^n 2 ^2 — 15 cos 2 + 10nin2 cost — nlj ■ 


(32) 

(33) 

(34) 


From the adjacent locations of several resonance centers in the i e plane (Figu re |T|) , it is natural to 
suspect that Chirikov’s resonance overlapping—a universal route to chaos (fChirikovj |1979 )—occurs in the 
mesh of the web. Stricto sensu, this was never rigorously demonstrated for the MEO problem, though [Ely| 
and Howell (1997) had derived all of the formulas necessary to do so. Appendix [A] gives the explicit formula 
for the lunisolar harmonic coefficients of the 29 resonant curves appearing in Figure [T] Recall that each 
center of the resonances (cf. Fig. [TJ are defined, in the i—e phase space, by the curves 

Cn = | (*, e) £ [0, 7^] x [0,1] : n ■ 0 = nig + n 2 h + = 0 j . 

The curves delimiting the maximal widths of each resonance are defined by 


W^ = <!(i,e)e[0,-]x[0,l] 


Note that the multiplication by d 2 T-L aec /dX 2 


0 = nig + n 2 h + n^T = ±2y^ 


d z H s 


dX 2 


x=x t 


(35) 


\x=x„ 


in Eq. 35 is due to the timescale operated in Eq. 25 


For each resonance, we compute the half-width excursion in the inclination and eccentricity phase space, 
according to Eq. |35| Figure [2] shows the centers and widths of the resonances in the region of semi-major 
axes between roughly 3 and 5 Earth radii, the range of validity of the present theory. At a = 19, 000 
km, the resonances are thin and the regions of overlap are quite narrow, indicating that regular orbits 
should dominate the phase space. As the ratio of the semi-major axis of the orbits of the satellite and the 
perturbing bodies (the perturbation parameter, e = e(a/aj\f)) is increased, the resonance centers for s j=- 0 
spread out and the libration regions become wider. We can observe a transition from the stability regime 
at a = 19,000 km to a Chirikov one where resonances overlap significantly at a = 29, 600 km. It is well 
known that mean-motion resonances become wider at larger eccentricity (q.v. [Mor bidelli 20021: the width 
of the lunisolar secular resonances exhibit a much more complicated dependence on the eccentricity, as well 
as the inclination. 

The Chirikov resonance-overlap criterion forms the extent of our theoretical analysis. This empirical 
criterion gives significant qualitative and quantitative predictions about the regions in action space for which 
chaotic orbits can be found; yet, it contains several limitations. It lacks a rigorous theoretical grounding 


and, as noted by Morbidelli and Guzzo (1996), it is not well tested on dynamical systems with many degrees 
of freedom. Chirikov’s geometrical argument neglects the coupling of the resonances—to say nothing about 
the deformation of their separatrices—and the role played by secondary resonances, and, as a result, the 
criterion often underestimates the threshold of transition from order to chaos. For a more complete and 
detailed analysis of the phase-space stability, we must turn to numerical explorations. As such, we use 
the fast Lyapunov indicator in the following section. As a systematic study of the entire parameter space 
represents a formidable task with significant computational requirements, we chose to focus on a reduced 
portion of the phase space that is of considerable practical interest for the navigation satellite constellations 
(see Figure [3|. 


The form of Eq. |34| contained in Ely and Howell (| 19971 is missing a factor of 2 in the second term in the brackets. 
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3 Numerical exploration of the phase space stability 

This section contains a numerical survey of the dynamical structures appearing in the MEO region by 
presenting an atlas of stability maps. Following a parametric approach, our main goals are (i) to give the 
geometry and extent of the stable, resonant, and chaotic domains, and (ii) to obtain a global picture of the 
resonant interactions on the emergence of chaos. In addition to quantifying the degree of hyperbolicity of 
generic orbits, we estimate the barriers of predictability by computing Lyapunov times. Furthermore, we 
demonstrate numerically that resonances and chaos are associated with transport in the i—e phase space. 


3.1 The numerical detection of the resonance overlapping regime 


In the past few decades, numerical investigations have played a wonderfully key role in studies on the 
long-term stability of dynamical systems. “The symbiosis between mathematical results and numerical 
computations,” writes Morbidelli and Guzzo (19961, “is very promising for the future developments of 
applied dynamical system science and, in particular, for Celestial Mechanics.” The fundamental work of 
Morbidelli, Guzzo, Froeschle, and others, have brought to light the existence of specific structures in the 
phase space when Nekhoroshev’s theorem is satisfied, which in turn imply this celebrated long-time stability 


result ( Morbidelli and GiorgilIi|1995 Morbidelli and Guzzo|1996| Morbidelli and Froeschle|1996 Morbidelli 

2002). Leaving aside mathematical intricacies, they essentially showed that even if chaos exists in the 


phase space (which is not precluded by the Nekhoroshev theorem), there always exists in a mesh of the 



Fig. 2: Lunisolar resonance centers (solid lines) and widths (transparent shapes) for increasing values of 
the satellite’s semi-major axis. This plot shows the regions of overlap between distinct resonant harmonics. 
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Fig. 3: Zoomed-in 


portion of Figure pi showing where 




we concentrate our numerical calculations. 


resonant web a no-resonant domain, which is filled by many invariant tori. Conversely, if the system is not 
in Nekhoroshev form, such no-resonant domains cannot be defined: resonances overlap and invariant tori 
become rare (Morbidelli and Giorgilli|l995 I. 

The divergent behavior of trajectories can be easily investigated numerically using the broad family of 
Lyapunov and affiliated indicators. The problem of the eventual overlapping of resonances is thus reduced 
to the numerical computation of indicators that distinguish between stable, resonant, and chaotic orbits. 
For such investigations, we chose from the chaos toolbox the fast Lyapunov indicator (FLI) (q.v. Froeschle 
et al. 1 997] ). Writing the n-dimensional dynamical system in first-order autonomous form, x = f(x), where 
x G R" and / : R" —I R” represents the vector field, the variational system in R 2n can be stated as 


x = f{x) 

■ fdf(x)\ 


(36) 


where w G R” stands for the deviation (or tangent) vector. Many first order stability indicators are based 
on the propagation of the variational system and on the monitoring of the stretching of w with time (see 
Skokos 2010 and references therein for a nice survey). The FLI, introduced in Froeschle et al. (1997|, 
follows from the variational system and enables the discrimination between ordered and chaotic motions 
(Froeschle and Lega||2000 Guzzo|[2002 1. The indicator at time t is defined by 


FLI(t) = sup log 11 w(t) 

T<t 


(37) 


where || • || denotes the Z/ 2 -norm. More rigorously, we should note FLI(t, xo, wo) instead of FLI(t) since 
the FLI computation requires the initialization of the variational system (Eq. 361. The choice of the initial 


tangent vector is a recurrent question when dealing with first order variational indicators, as discussed by 
Barrio et al. (2009). The FLI atlas—a collection of FLI maps—presented here depend on two parameters: (i) 


the initial semi-major axis, ranging from roughly 3 to 5 Earth radii (the prominent MEO region) and (ii) a 
hierarchy of physical models, in order to illustrate the different dynamical mechanisms and how resonances 
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Table 1: Perturbations added to the central part of the geopotential for the numerical stability analysis. 
We refer to Appendix [C] for more details. 


Dynamical models 



Zonal terms 

Tesseral terms 

Lunar perturbation 

Solar perturbation 

model 1 


J 2 

not considered 

up to degree 2 

not considered 

model 2 


j 2 

not considered 

up to degree 2 

up to degree 2 

model 3 

J 2 , 

5 J5 

not considered 

up to degree 4 

up to degree 3 

model 4 

j 2 , 

,J5 

up to degree and order 5 

up to degree 2 

up to degree 2 


and chaos appear. The physical models considered are schematically summarized by Table [l] and discussed 
in more detail in Appendix [C] along with a description of the numerical setup. 


3.2 FLI stability atlas: hyperbolicity and predictability 

Figures [4] through [7] shows the results of the FLI analysis corresponding to the zoomed-in portion of the 
action space presented in Fig. [3] for the four dynamical models of Table [l] On each of the pictures, the 
initial conditions io, eo are associated to the corresponding FLI value using a black-blue-red-yellow color 
scale. The FLI of all regular orbits have approximately the same value and appear with the same dark 
blue color. Light blue corresponds to invariant tori, red and yellow to chaotic regions, and white to collision 
orbits. The striking feature of Figs.[4}j7]is the transition from ordered to chaotic motion as the perturbation 
parameter increases. For ao = 19,000 km, ordered motions dominate in the phase space. Resonant structures 
are mostly isolated or emanate in a “V” shape at the birth of inclination-dependently-only resonances at 
zero eccentricity. Chaotic orbits are mostly confined along the principal V shape of the critical inclination 
resonance (63.4°) and to small pockets of instability at higher eccentricities. As the semi-major axis increases 
to 24, 000 km, the geometrical organization and coexistence of regular and chaotic orbits becomes more 
complex. At 25, 500 km, more and more of the resonant tori are rendered unstable, and the width of the 
stochastic layers and their degree of overlap are increased. A common characteristic between these maps 
and the 24, 000 km case is the location of re-entry orbits along the 56.1° resonance. Finally, at ao = 29, 600 
km (the nominal semi-major axis of the Galileo constellation), roughly speaking, all orbits belonging to the 
domain delimited by (io,eo) £ [52°,71°] x [0,0,5] are resonant or chaotic. Re-entry orbits display a very 
complicated geometry and may now concern quasi-circular orbits near i o = 55°; that is, strong instability 
exist that can potentially affect Galileo-like orbits. 

A fundamental conclusion reached via this parametric approach, using a hierarchy of dynamical models, 
is that model 2 can be legitimately declared as the basic force model; i.e., the simplest physical model 
that can capture nearly all of the qualitative and quantitate features (dynamical structures in the maps, 
degree of hyperbolicity, domains of collision orbits, etc.) of more complicated and realistic force models. 
In particular, there are no significant changes in the FLI maps of Figs. [6] and [7] when compared to that 
of Fig. [5] Force model 2 differs from model 1 only by the presence of the solar third-body perturbation, 
developed at the quadrupole order. Comparing Figs. [5] and [4] it is clear that that solar perturbations play a 
non-negligible role on the long-term dynamics. We note, specifically, the manifest widening of the resonant 
regions and the increase of chaotic orbits near the inclination-dependent-only resonances, principally near 
■i = 56.1° and 63.4° for all eccentricity values. Moreover, the volume of collision (re-entry) orbits is larger 
when solar perturbation are taken into account, as is well illustrated in panels (b), (c), and (d) for highly 
eccentricity orbits along i = 69°. This numerical finding confirms, a posteriori, the analytical refinement 
given by Eq. [24] The robustness of the physical model 2 was further tested at a smaller phase-space scale, 
as shown in Fig. [8] giving a magnification of a zone containing many secondary and transverse structures 
for ao = 24,000 km. To show finer details and sharper contrast, the resolution and propagation time were 
enhanced for an objective comparison between the various force models. We can see that, even at this 
scale, force model 2 can unquestionably be considered the basic physical model: changes in the dynamical 
structures, and the like, with models 3 and 4 are nearly undetectable. These maps (Fig. [8] also reveal the 
extraordinary richness of the phase space and the apparition of transverse and thin hyperbolic manifolds 
(structures with high FLI values). It is worth nothing that re-entry orbits, appearing when the solar 
contribution is included, are located precisely along those thin manifolds, already present with force model 
1. Furthermore, the main collision region in the vicinity of the 56° inclination-dependent-only resonance 
widens significantly with model 2. The identification of the relevant dynamical model is undoubtedly the 
first question to be addressed when dealing with real (physical) problems. Strangely enough, the isolation of 
the basic force model for the MEO problem, was largely missing from the literature, though was speculated 
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on in Rosengren et al. (2015). Thus, all future work in this area can be performed under this basic model 
of oblateness and lunisolar perturbations, and, without loss of generality, the obtained results may be 
considered representative of the other more refined, and complicated, models. 

It should be emphasized that all of these charts have been obtained by varying only the initial inclination 
and eccentricity, with no regard for the initial phases, themselves being fixed for all computed FLI. Figures[9] 
and [To] illustrate the inherent difficulty to capture the dynamics of a six-dimensional phase space in a plane 
of dimension two (recall that Hamiltonian, Eq. 14 is three DOF) (Ric hter et al.||2014 |. They present the 
FLI maps for ao = 29, 600 km under model 2, at two different phase-space scales, in which the results of the 
latter may be considered representative of Galileo-like orbits. We can observe how the dynamical structures 
(stable, resonant, chaotic, or collision orbits) evolve by changing the initial angles cu. 12, or even the initial 
configuration of the Earth-Moon system (equivalent to changing the initial epoch of the simulation). Despite 
this dependence on the initial phases, however, the regime of the global phase space still persists: the region 
is dominated by chaotic orbits. This does not appear to be the case for the zoomed-in portion of the phase 
space, where it is clear from Figs. 10(a) and 10(d) how the initial epoch may strongly influence the stability 
analysis. Certainly, the dependence of the dynamical structures on the initial phases (u;, 12, 12 m) is a subject 
that requires further studies, and the associated difficulties are, in essence, due to the representation of the 
dynamics in a reduced dimensional phase space (Richter et al.|[2014|. 




50 55 60 65 70 75 

i [deg] 


(a) ao = 19, 000 km. 


(b) ao = 24, 000 km. 



(c) ao = 25, 500 km. 


(d) ao = 29, 600 km. 


Fig. 4: FLI stability maps for dynamical model 1. 


Besides the quantification of the local hyperbolicity, the predictability is an another important aspect of 
dynamical interest that we would like to discuss here. The Lyapunov time T^, defined as as the inverse of 
the maximal Lyapunov characteristic exponent (mLCE), represents physically the barrier of predictability 
in the dynamical system. Recall that this time, defined by 


T r = 


lim 

t—¥ OO 


log||w(r)|| 


(38) 


goes to infinity for regular orbits (because the mLCE decreases to zero), and, as a result, the dynamical 
system is predictable for all future time t. However, for chaotic systems, the mLCE converges to a positive 
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(c) ag = 25, 500 km. 


(d) ao = 29, 600 km. 


Fig. 5: FLI stability maps for dynamical model 2. 


value and converges to a finite value. We have estimated Eq. [38] numerically for the charts of Fig. [lO] 
and we have found that many of the orbits have Lyapunov times on the order of decades, as shown by 
Fig. [ll] for a specific case. For the propagation time chosen, slightly larger than 530 years, the most stable 
orbits in the maps exhibit, roughly speaking, a Lyapunov time of between 100 and 120 years. Thus, in 
Fig- HU the orbits with a Lyapunov time bigger than 100 years have been assigned to an orbit with a 
Lyapunov time equals to 100 years. This does not alter the stability results, but only provide more visual 
contrast in the maps. This order of predictability, found to be very small indeed for Galileo-like orbits, 
is a nice invitation to the general reflection of the practical meaning of long-time simulations, so often 
preformed in astrodynamics. In fact, the problem should be considered as probabilistic in nature and the 
evolutions of orbits, far beyond the Lyapunov time, have to be considered as only a possible future. In this 
sense, the long time tracking of individual orbits appears to be inappropriate. The problem requires more 
fundamental statistical studies that describe the properties of an ensemble of trajectories: this approach is 
far from common in our community. 


3.3 Transport in phase space 


The local hyperbolicity and predictability in the MEO region, both first-order variational quantities, have 
been quantified via FLIs and numerical estimation of Lyapunov times, respectively. Another important 


feature associated with chaos, which has yet to be addressed, is the transport properties (Todorovic et al. 


2008); that is, the possibility for the dynamical system to explore its phase-space domain. As typified by 


the asteroid Helga (Milani 19921, stable chaos may occur in actual (physical) solar system dynamics, for 
which orbits with short Lyapunov times display no signifiant changes on very long timescales (even several 
orders of magnitude larger than T^). For the MEO problem, stable chaos is evidently not at the heart 
of the instabilities and, on the countrary, the apparent growth of the eccentricity on short timescales was 


how chaos was initially suspected (Rosengren et alJ2015). The most rewarding suggestions of the transport 
properties concerning MEO dynamics are due to Ely["(2002), who showed the tendency for typical GPS 


orbits to follow the resonant skeleton give by Fig. [l] This explains how quasi-circular orbits may become 
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(c) ag = 25, 500 km. 


(d) ao = 29, 600 km. 


Fig. 6: FLI stability maps for dynamical model 3. 


nearly hyperbolic on millennia timescales. Similar results were obtained recently by Rosengren et al. (2015 I, 


using stroboscopic techniques, who demonstrated this phenomenon on much short timescales for the various 
navigation constellations. Here, we use our FLI analysis, which more clearly reveals the extent and geometry 
of the chaotic sea in the phase space, to further enliven this idea, by showing numerically how the dynamical 
structures in the FLI maps affect the long-term evolutionary properties of the motion. 


Figure 12 shows the long-term evolution (plotted stroboscopically once per lunar nodal period) of 


orbits with initial inclinations and eccentricities in particular regions of the phase space superimposed on 
the background dynamical structures obtained from the FLI computations. For orbits initial located far 
from the resonances, nothing dynamically interesting happens in terms of transport: the motion is quasi- 
periodic and the excursions in the action space are modest, bounded, and confined by the invariant tori. 
Figures |l2(a) 12(c)| show that in the vicinity of isolated resonances (where no adjacent resonances exist), 
the orbits have the tendency to explore a larger portion of their phase space, evolving along this resonance. 
The evolutionary properties of the eccentricity and inclination are mainly determined by the shape of the 
resonance, as revealed by the FLI analysis. Because of the geometry of the dominate isolated resonances 
(the V shapes clearly distinguishable for ao = 19,000 km), we see how some orbits can have a very stable 
inclination despite significant growth in eccentricity (evolving along the resonant line). For orbits with 
initial conditions in the resonance overlapping regime, where invariant tori are destroyed, we find that 


excursions of the eccentricity and inclination are fast and macroscopic (Fig. 12(d)). The exploration of a 
large phase-space volume is permitted, and, despite the erratic appearance of the motion, when sampled 
stroboscopically we see a tendency for the orbits to jump from one resonance domain to another, being 
confined to the chaotic sea. Figure fl2(d)| also shows the evolution of orbits in the stable domains (marked 
green and yellow), where the variations are quasi-periodic. Such results were further confirmed by a large 
amount of simulations, and thus we can conclude that the FLI maps clearly reveal how transport is mediated 
in the phase space. 
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(c) a o = 25, 500 km. 


(d) ao = 29, 600 km. 


Fig. 7: FLI stability maps for dynamical model 4. 


4 Testing the Chirikov criterion 

The dynamical structure of the phase space of the MEO region was enlightened using two different, but 
complementary, approaches: the heuristic Chirikov criterion and a CPU-expansive numerical stability sur¬ 
vey. We would like to stress at the outset that, to the best of our knowledge, the Chirikov principle has 
never been applied to a real (physical) system of more than 2 degrees of freedom; recall that in our study 
the autonomous version of the Hamiltonian is 3 DOF, leading to a 6 dimensional phase space, whereby the 
resonance widths involve two actions. The realistic nature of Chirikov’s analytical results is then of primary 
importance, and so a synthesis of the theoretical and numerical aspects of this work is presented here. 

Figure [13] presents the lunisolar secular resonance centers and widths of Fig. |3j now appearing as light 
gray, superimposed on the corresponding FLI stability maps of Fig. [5j both obtained under the same basic 
force model of oblateness and lunisolar perturbations. For a o = 19,000 km, where the resonances are mostly 
isolated (the fundamental assumption made in our analytical treatment), the agreement is both qualitatively 
and quantitatively very satisfactory. The FLI numerically detected V shapes, emanating from 56.1°, 60°, 
and 63.4° at zero eccentricity, are well captured analytically. As the semi-major axis increases, however, the 
discrepancy with the isolated resonance hypothesis becomes manifest, as was noted from the FLI analysis. 
Accordingly, the comparison between the two approaches becomes more and more difficult to quantify, 
though we note that the most important qualitative feature—the fact that the systems tends towards the 
overlapping chaotic regime—can be seen from both. The analytical development captures indubitably the 
tendency of global overlap, even if a 1 — 1 correspondence between the overlaps and the numerical detection 
of chaos is hard to formulate. Nevertheless, we note that some numerical structures are nicely captured 
analytically, such as the widths of the resonances for moderately eccentric orbits near the inclination of 
54°-60° and 62°-64° for ag = 24, 000 km, or near the critical inclination for a o = 25, 500 km. Also apparent 
for all semi-major axes is the large number of chaotic or escaping orbits along the 69° resonance at high 
eccentricities at the intersection of multiple resonant curves. 

Reasons for any discrepancies surely arise from the fact that the Chirikov criterion is imperfect by 
nature: resonances are treated solely in isolation. Certainly, more sophisticated and rigorous analytical 
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(c) dynamical model 3 


(d) dynamical model 4 


Fig. 8: Zoomed-in portion for ao = 24, 000 km near the 2cL> + 12 only-dependent-inclination resonance under 
the various dynamical models. Initial conditions have been propagated from the initial epoch 2 March 1969 
until the final date set to 15 November 2598. The precise detection of the stable manifolds allows to predict 
the set of re-entry orbits. These maps also corroborate model 2 as the basic physical model. 


methods that treat resonance interactions, such as those detailed in Haller (1999), should be pursued in the 
future. But despite the approximations of the Chirikov criterion, it has permitted us to gain a formidable 
intuition into the nature of this physical problem. The application of this analytical criterion is a strong 
testimony that despite the great power and scope of modern computers, there is still a place in celestial 
mechanics for pen-and-paper calculations in the style of our forebears. 


5 Discussion and future work 


We present below several interesting aspects of the MEO problem that have not been adequately addressed 
here; points on which to concentrate in future works. 

1. We have pointed out the critical role of the initial phases (w, 12, 12 m) on the dynamical structures in the 
FLI maps, leading to different chaotic pictures in the phase space. Yet, we gave no information about 
which angles will ensure or avoid chaos for a given initial inclination and eccentricity. This question is 
clearly of significant practical interest and needs further investigation. 

Even if we demonstrated a 1-1 correspondence between chaos and large-volume exploration of phase 
space, the precise nature of the transport has not been specified. In particular, no distribution of the 
increments Ae with time was given. Such analytical modeling of f(Ae, t ) is currently being studied, 
based on considerations in Murray et al. (1997) or Varvolglis (2004). 


2 . 


3. The general link between local hyperbolicity (chaos) and stochasticity (Chirikov 1979) needs to be 
studied, a question which has hitherto remained seriously underrated in celestial mechanics. Statistical 
descriptions of the action evolution, generally through adequate averaged quantitates, are at the foun¬ 
dations of statistical mechanics and ergodic theory—the theory of the long-term statistical behavior 
of dynamical systems. The injection of the concepts and tools of these fields to the applied (physical) 
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Fig. 9: FLI stability maps for ao set to 29,600 km under force model 2. Unless otherwise stated, the initial 
epoch of the simulations, which determines the initial dynamical configuration of the Earth-Moon-Sun 
system, is 02 March 1969. 


MEO problem must still be performed. This need, moreover, is deeply stimulated by the estimated short 
Lyapunov times, implying that predictions must be statistical in nature. 

4. In the Chirikov overlapping regime (ao = 29,600 km), re-entry orbits represented a significant portion 
of the phase space. The exact period of time after which an orbit collides, however, is not illustrated 
in the FLI maps, so that there might exist orbits that collide after 100 years or after 500 years, both 
represented in the maps by the white color. We preferred here to keep our study succinct and focused 
only on dynamical principles, but such considerations are obviously relevant for the analysis of disposal 
strategies for the four constellations located in this precarious region of space. 

5. Though averaging methods have proven their use in celestial mechanics, given the criticism of Arnold 
to the so-called “averaging principle” that was in vogue 50 years ago, we are currently testing the 
simulation of the whole system (the non-averaged system) to test the validity of the averaging procedure 
with respect to chaotic dynamics (Arnold'et al.[~2006; Daquinp0l51 


6 Conclusions 


We have presented a 2.5-DOF secular Hamiltonian governing the long-term dynamics of the MEO region 
under oblateness and lunisolar perturbations. Using this analytically tractable Hamiltonian, written in 
terms of the Delaunay elements and reduced to a 1-DOF pendulum system near specified resonances, we 
have estimated the widths of the dense network of secular lunisolar resonances permeating the i—e phase 
space. As the semi-major axis recedes from Earth, scanning the MEO region (from 3 to 5 Earth radii), 
we found a transition from a globally stable regime, where resonances are thin and well separated, to a 
globally unstable one, where resonances widen and overlap. Surprisingly, the application of the well-known 
analytical Chirikov criterion, frequently used for architectural celestial mechanical studies, was never fully 
applied to the MEO problem, because the fundamental work of Ely and Howell (1997]) was overlooked 
by the community for nearly twenty years. This analytical methodology was completed here, and tested 
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Fig. 10: Zoomed-in portion of Fig. |9j representative of Galileo-like 
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Fig. 11: Map of Lyapunov-time expressed in years. The Lyapunov time represents the barrier of predictabil¬ 
ity of the dynamical system. Far beyond this time, the system loses the memory of his initial conditions 
and statistical properties are more suitable to characterize the evolutionary properties of the motion. Initial 


conditions are those of Fig. 10(a) 


against predictions from a semi-numerical stability analysis, which focused on the region of MEO populated 
by the navigation constellations. In essence, both approaches agree and confirm the transition to chaos; 
yet it is well known that the approximation of treating the resonances in isolation breaks down in the 
overlapping regime, implying that the quantitative extent and geometry of the chaotic domains can only 
be obtained numerically. Our parametric and extensive numerical simulations, moreover, have permitted 
the isolation of the basic physical model governing the dynamics in this region. Specifically, we have shown 
that the grandiloquent force models, often employed in such studies, are not required to capture the 
stable or chaotic features of the phase space. Future analytical or numerical effects can be made with 
the basic 2.5-DOF averaged model of oblateness and lunisolar perturbations, developed to quadrupole 
order. The semi-analytical investigations have also illustrated the extraordinary richness of the dynamical 
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(b) transport along the isolated 2lj + ft = 0 resonance; 
zoomed-in portion. 



(c) Transport along the isolated critical inclination resonance. 
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Fig. 12: Illustration of the transport properties for orbits with initial condition near isolated resonances or 
for orbits with initial condition in the resonance overlapping regime or far from them. All orbits are semi- 
analytically propagated using the basic force model 2 and are stroboscopically overlaid onto the FLI maps 
obtained with the same force model. The nature of the resonances compel the transport in the phase space. 
The structures in the FLI maps, in addition to quantifying the hyperbolicity, also reveal the preferential 
and possible routes for transport in phase space. 


structures in the i—e phase space, whose full understanding is far from complete. The predictably of typical 
navigation satellites has also been quantified through the estimation of Lyapunov times, finding a barrier 
of predictability of only decades. It is hoped that such results will stimulate and guide the community 
towards systematic and fundamental statistical approaches, more suitable for describing the distribution of 
the actions. We have also investigated the transport characteristics in the phase space, and how resonances 
and chaos influence the evolutionary behavior of the inclination and eccentricity. The nature of this transport 
will be discussed in a forthcoming paper, as with some general stochastic properties of the motion induced by 
the hyperbolic (chaotic) nature of the dynamical system. The effective application of these results towards 
the management of the navigation satellite systems is a challenging task, the goal of which is to exploit 
the instabilities to identify suitable collision orbits or to use the associate transport routes to reach stable 
regions in phase space (parking orbits). In this regard, we stress that analytical and numerical techniques 
cannot be separated, and a complete, logically ordered picture can be obtained only by the application of 
both methods jointly. 
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Fig. 13: FLI stability maps versus Chirikov’s resonance overlap for force model 2. 
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A Special functions in the lunisolar disturbing potential expansions 

Recall that the lunar disturbing function can be written in the compact form 

2 2 2 

TZm = EEE ^2—2 p,m,±s COS < &2—2 p,m,±s? 

m =0 s=0 p=0 

with harmonic angle and associated harmonic coefficient 
®2-2p,m,±s = (2 - 2 p)u + mQ ± s(J?M tt/ 2) - y s 7T, 


, M 
ho- 


rKvTH (-l)L m / 2 Je T n#-^^, m , P WF 2 . S 4(*M)^2 > p, 2p -2(e)(-l) m+ ^r TS (e), 


2-2 P ,m,±s a ^ (1 _ e ^ )3/2 v -y " m (2 + m)\ 

so that a lunar resonance occurs when 

ip2-2 P ,m,±s = (2 - 2p)cl) + mf2 ± sfi M ~ o. 
The solar disturbing function can be written in the compact form 

2 2 

K s = EE h2-2p,m COS <?2-2p ,ra i 

m =0 p=0 
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with harmonic angle and coefficient 

^ 2 - 2 p,m = (2 - 2 p)uj + m(n - l?s), 


hL 


2 p,m 


MS a 


x|(l — e l) 3/2 (2 + m)! 


(2 — m)\ 

SmTTl - -(r-F2,m,p(*)-F2,m,l(*s)-f^2,p,2p-2(e). 


A solar commensurability occurs when (2 — 2p)u + ml? « 0, or equivalently when ip 2 - 2 p,m,o ~ 0. We give 
here the the explicit formula needed to calculate the widths of the 29 distinct curves of secular resonances, 
six of which are locations of lunisolar resonance (i.e., the inclination-dependent-only cases), appearing in 
Figure [l] 

The general Hansen coefficient x!’ m (e), which permits the full disturbing function (prior to averaging) 
to be developed in terms of the mean anomalies, is a function of the orbit eccentricity and is given by the 
integral ( Hughes| |l980) 


-i /*27T l 

xj: m (e) = — I 0 cos {mf - kM) d M. 


(39) 


For k = 0, exact analytical expressions exist for the zero-order Hansen coefficients X l Q ’ m for all values of l 
and m. The integral for X l 0 m can be written as 


/ 1 r 271 /r\ l 

A » -^,0 


/•27T l 

J cosm/dM. 


(40) 


Note that since cosine is an even function, it is only necessary to obtain expressions for m > 0 as X l Q ’ m = 
Xg~ m . For l > 1 and 0 < m < l, the integrals can be evaluated as 


V -Z,m ( e.\ nl f l -\- m ^ (m l 1 m l 2 i 

A o = V 2/ { rn - 2 - + >’ 


(41) 


where F(q, /3, 7 ; x) is a hypergeometric function in e 2 . Several formulae of recurrence have been derived that 
greatly facilitate the calculation of these coefficients, however, for our purposes, H 2 lP , 2 p- 2 (e) = A'g ' 2 2p ( e ), 
for p = 0 , 1 , 2 , can be easily evaluated as 


t- 2,±2 b 2 xr2,0 -,.32 

A o = 2 e ’ A o = 1 + 2 e ' 

The Kaula inclination functions F 2 t m,p(i) and F 2 jSi i(*m) are given by (Kaula 1966) 

(4-2 1 )! 


^2 ,m,p(i) — ^ 
t 

m 

E 

s=0 


(2 — t)\(2 — m — 2t)\2 4 ~ 2t 


• 2—m—2t • 

sin 1 


*E 


2 — m — 2t-\-s \ 1 m — s 
p — t — c 


(-1) 


c—k 


(42) 


(43) 


where k is the integer part of (2 — m)/ 2 , t is summed from 0 to the lesser of p or k, and c is summed over 
all values making the binomial coefficients nonzero. Expressions for F 2 , m ,p(i) up to m = 2, p = 2, are given 
in Table [2j 

The expression for the Giacaglia functions, computed from Eq. [4j are listed in Table [3] 

The commensurability condition, harmonic angle, and associated harmonic coefficient for all resonance 
center curves are given in Tables [4] and [5] color coded to match that of Fig. [l] Note that when the angular 
arguments are the same (excepting a constant phase), as in each of the five resonance curves stemming 
from the critical inclination (63.4°) at e = 1 or the lunisolar inclination-dependent-only resonances, the 
associated harmonic coefficients must be combined into one single term for the calculation of the resonant 
width, as detailed in Section [2.2| 
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Table 2: Kaula’s inclination functions 


F 2 ,m,p{i) from Kaula (19661. 


m 

p 

F 2 ,m,p(i) 

0 

0 

—3 sin 2 i/S 

0 

1 

3 sin 2 i/4 - 1/2 

0 

2 

—3 sin 2 i/S 

1 

0 

3sinz(l + cosz)/4 

1 

1 

—3 sin i cos i/2 

1 

2 

— 3sinz(l — cosz)/4 

2 

0 

3(1 + cosi) 2 /4 

2 

1 

3 sin 2 i/2 

2 

2 

3(1 — cosi) 2 /4 


Table 3: Giacaglia’s tesseral harmonics rotation functions [/™ ,=FS for the Moon, where C = cos(e/2) and 
S = sin(e/2). 


m 

0 

-1 

1 

-2 

2 

0 

1 - 6C 2 + 6C 4 

-2CS" 1 (2C 4 - 3C 2 + 1) 

—2CS(1 - 2C 2 ) 

C 2 S~ 2 (C 2 - l) 2 

C 2 S 2 

1 

—3CS _1 (2C 4 — 3C 2 + 1) 

S -2 (4C 6 - 9C 4 + 6C 2 - 1) 

C 2 (4C 2 - 3) 

-CS~ Z (C 2 - l) 3 

-C 3 S 

2 

6 C 2 S~ 2 (C 2 - l) 2 

—4CS i- 3 (C' 2 - l) 3 

—4C 3 S _1 (C' 2 - 1) 

S~ 4 (C 2 - l) 4 

c 4 


B Inclination-dependent-only lunisolar resonances 


Figure 14 shows the lunisolar inclination-dependent-only resonance centers (solid lines) and widths (trans¬ 
parent lobes) for increasing values of the satellites’s semi-major axis, computed with the solar contribution 
and without (shaded gray and without an edge color). The solar harmonics, neglected in previous works 
(Ely and Howell 19971, play an increasingly important role as the semi-major axis is increased, either by 
widening or shrinking the resonance domain. The lobe shape of the lunisolar apsidal resonances (q.v. B reiterj 
2001) and the rocket shape of the nodal resonance are easily explained by the dependence of the widths on 
the orbit eccentricity. The widths depend on e through their associated harmonic coefficients (Appendix |A|, 
which comes in through the zero-order Iiansen coefficients, and through Eq. [33] From Eqs. |26[ [3] [TJ and 
42] v —> 0 as e —> 0 for p = 0 or p = 2, but not for p = 1. To extend the result to the case where e —» 1, let 
us note that the equation 


mw + n 2 fi = ±2^v.d 2 H sec /dX 2 \ x=x * 


(44) 


is equivalent to 


nP a (i) = ±2u(e)Jv.d 2 'H sec /dX 2 \ x=x ^ 


(45) 


where k a constant (independent of e and i), u(e) = (1 — e 2 ) 2 and P n (i) = ni(5cos 2 i — l)/2 — n 2 cos i. When 
e —> 1, the right-hand side of Eq. (451 goes to zero, and, as a result, i must be a root of P n in the limit 
e —> 1. 


C Numerical setups 

The dynamical model we used in the propagation of the orbits accounts for the perturbations stemming 
from the Earth, the Moon, and the Sun. To simplify the mathematical problem, all of these gravitational 
perturbations have been analytically averaged over the mean anomaly M of the satellite (the fast variable), 
and propagated using numerical integrations, which is much faster thanks to the absence of short-periodic 
variations. This is a well known and efficient semi-analytical approach to study the qualitative evolution 
of orbits over very long timespans. The averaging procedure has been done for the Earth’s geopotential up 
to J5. For tesseral resonances (located for certain semi-major axes), where there exists a commensurability 
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Table 4: Lunar secular resonance conditions, harmonic angles, and associated harmonic coefficients needed 
for the resonant width computations of Section 


2.2 


Lunar commensurabilities 


■02 — 2 p,m,±s 

(f)M 

2 — 2p,m,zbs 

/,M 

fc 2 —2p,m,is 

02,2,0 

2oj + 2ft 

15mm“ 2 (3s 2 i m - 2)C 2 5“ 2 (C 2 - l) 2 2 2 

324(1 -*)■/■ ' (1+a) 

02,2,1 

2cj -(- -|- iT2]yj — 7T 

15/lM“ 2 S*MCiMC'5 _3 (C 2 - l) 3 2 2 

.„3,i 2 13/2 e C 1 +“) 

!6<(1 - eh) 3 ' 2 

0 2,2,-1 

2cj -|- 2i2 — 

15/lM« 2 S*MCiMC ,3 S _1 (C 2 - 1) 2 2 

.,.3,1 2 13/2 e (! +“) 

!6<(1 - eh) 3 ' 2 

02,2,2 

2CJ 2i2 2i2]y[ — 7T 

l5n M a 2 s 2 i M S~ /L (C 2 - l) 4 | 

64^(1 “ e M) 3/2 

■02,2,-2 

2cj -|- 2i2 — 2/2jy[ 7r 

15MM« 2 s 2 i M C 4 , n . -i2 

Ma 3 M (l-e 2 M ) 3 ' 2 ( + C) 

■02,1,0 

2oj + ft 

15/r M a2 (3s 2 i M — 2)C5 _1 (2C 4 — 3C 2 -(-1) 2 , 

3 /, 2 13/2 e SH 1 +«) 

ISa^l - eh) 3 ' 2 

■02,1,1 

2cj -|- .T2 if2]yi — 7T 

15/iM« 2 s*MciMiS ,_2 (4C 6 — 9C 4 + 6C 2 — 1) 2 .. 

1C 3 2 13/2 e s *( 1 + c *) 

- e M) 3/2 

02,1,-1 

2 Cel “I - i(2 — 

15/i M a 2 siMciMC' 2 (4C 2 - 3) 2 Yl . .. 

ic 3(1 2 13/2 e s *( 1 + «) 

- e M) 3/2 

02,1,2 

2cj -|- .T2 -|- 2 — 7r 

15// M a 2 s 2 i M C,S _3 (C 2 - l) 3 2 Y1 

-- - e sz(l + cz) 

16“m(1 - eM) 3/2 

02,1,-2 

2CJ -|- — 2 -|- 7T 

IS/im^s 2 ^^ 3 ^ 2 . .. 

-s—-„ NO e si 1 + ci) 

160^(1 - eh) 3 ' 2 

02,0,0 

2cj 

15// M a 2 (3s 2 i M - 2)(1 - 6C 2 + 6C 4 ) 2 2 ■ 

64<(1 - e 2 M ) 3 ' 2 

02,0,1 

2cj iT2]y[ — 7r 

45fiMe 2 siMciMCS~ 1 (2C 4 — 3C 2 + 1) 2 2 
32<(l-e 2 I ) 3 / 2 

02,0,-1 

2cj — i/ M 

45//M a2 s*MCi]y[C>S : (l — 2C 2 ) 2 2 

32^(1 -e 2 ,) 3 / 2 eS? 

02,0,2 

2cj + 2J?]VI — 7T 

45/tM« 2 s 2 *MC' 2 S _2 (C 2 - l) 2 2 2- 
64*4(1 - e 2 M ) 3 ' 2 

02,0,-2 

2cj — 2i2]yi + 7r 

45fi M e 2 s 2 i M C 2 S 2 22 . 

64 a M(l — e M) 3 / 2 * 

0-2,0,0 

—2uj 

15// M a 2 (3s 2 i M - 2)(1 - 6C 2 + 6C 4 ) 2 2 ■ 

64a 3 I (l - e 2 ,) 3 / 2 

0-2,0,1 

—2lj -)- — 7 r 

45//Mffl 2 siMC/MC'‘S' _1 (2C 4 — 3C 2 + 1) 2 3 

Wahil-eh) 3 ' 2 

0-2,0,-1 

—2lj — 

45//Ma 2 siMC*MC'S(l — 2C 2 ) 2 2 
32<( 1 -e 2 1 ) 3 / 2 682 

0-2,0,2 

— 2biJ -)- 2/2fy[ — 7T 

45/t M a 2 s 2 j M C 2 S _2 (C 2 - l) 2 2 2 - 
64<(1 - e 2 ,) 3 / 2 

0-2,0,-2 

— 2(jJ — 2/2]y[ -(- 7T 

45/iM® 2 s 2 *MC 2 'S’ 2 2 2 • 
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Lunar commensurabilities (cont.) 


'4^2 — 2p,m,±s 

rZ>M 

^2 — 2 p,m,is 

lM 

L 2 — 2 p,m,is 

lp-2,1,0 

—2u + (2 

15/i m ® 2 (3s 2 i M — 2)CS — 1 (2C ' 4 — 3C 2 + 1) , 

3 2 M/2 esl1 “) 

i>-2,l,l 

— 2u) -j - £2 ~\~ — 7T 
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15fi M a 2 siMciMC 2 (4C 2 - 3) 2 

-t,- 77 ———-e sz( 1 — cz) 

16<(1 -e ^) 3 / 2 

02 , 1,2 
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15mm® 2 s 2 *mCS - 3 (C 2 - l ) 3 2 
-o-o- 7C77. -e sz( 1 — cz) 

16 “m (1 - e h) 3/ 

Ip- 2 , 1,-2 

— 2uj £2 — 2 £2^[ H - 7 r 

15//m« 2 s 2 *mC ,3 S 2 
-b- a— esz( 1 — cz) 

16«m( 1 — e M ) 3 / 2 

02 , 2,0 

— 2 oj -|- 2(2 

15/i M a 2 (3s 2 i M - 2 )C 2 S~ 2 (C 2 - l ) 2 2 2 

'“-“I 

02 , 2,1 

— 2uj -(- 2 £2 £2m — 7T 

15flMd 2 SImCImCS~ 3 (C 2 - l ) 3 2 2 

1 &>L <1 - 4)»/> 

02 , 2,-1 

— 2uj -(- 2 £2 — 

15MM“ 2 S*mCImC ,3 S ,_1 (C' 2 - 1 ) 2 2 

1«4(1 - ^ e< 1 

02 , 2,2 

— 2lj -\- 2 £2 -\- 2 £2^/i it 

15fi M a 2 s 2 i M S~ 4 (C 2 - l ) 4 2/1 -\2 

®4 q m( 1 — e M ) 3 / 2 e( 

02 , 2,-2 

— 2lj + 2£2 — 2 £2^/i ~\~ it 

15^ M a 2 s 2 i M C 4 2 2 

64 “m ( 1 ~ e u) * 

" 00 , 1,0 

(2 

3// M “ 2 (3s 2 *m — 2)CS~ 1 (2C 4 — 3C 2 + 1) , , 2 , . . 

„ 3 ,, 2 / 3/2 (2 + 3e )sici 

0 >, 1,-1 

£2 — £2m 

3// M a 2 s*MciMC 2 (4C 2 - 3) 2 . . 

Q 3 /, 2 \3/2 2 + 3e SlC * 

^mC 1 - e h) 3/2 

0 ),l ,-2 

£2 — 2 £2yi -\- TT 

3 m aVi M C 3 S 2) ; 

8a3 I (l-e 2 I ) 3 / 2(2 + 3e )SJ " 

V 1 0 , 2,0 

2(2 

3/i M ® 2 (3s 2 JM - 2 )C 2 S~ 2 (C 2 — l ) 2 , 0 , Q 2 \ 2 ■ 

16 a M ( l — e M } 3 / 2 ^ ^ 

Ip 0 , 2 ,-l 

2 £2 — £2y[ 

‘i^M<P 2 si M ciMC :i S- 1 {C 2 -l) 2 / 2 - 

84(1-4)^ 12 + 3,1 

lp0,2,-2 

2£2 — 2£2m + TT 

3#iMH 2 8 2 i M C 4 2,2.. 

32<(l-e 2 I ) 3 / 2(2 + 3e )Sl 


between the frequency of the satellite’s mean motion and the sidereal time, a partial averaging method 
was applied to retain only the long-periodic perturbations. This is equivalent to retaining in the Earth’s 
geopotential only the slowly varying quantities associated with the critical resonant angle, as it is technically 
detailed in Morand (2013). Tesseral resonances are analytically well modeled by a 1.5-DOF Hamiltonian 
and they primary affect the semi-major axis of the satellite, leading to a thin chaotic response of the 
system (the 1.5-DOF bounded chaos is physically equivalent to intermittency phenomenon on the semi¬ 
major axis). Despite the fact that chaotic response of the system is very confined in phase-space (on the 
order of kilometers), we have included these perturbations in the full system to check if such effects may 
couple with the lunisolar effects on longer timescales. The coefficients of the Earth’s gravity field come from 
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Table 5: Solar secular resonance conditions, harmonic angles, and associated harmonic coefficients needed 
for the resonant width computations of Section 


2.2 


Solar commensurabilities 


1p2 — 2p,m 

2 — 2p,m 

h s 

n 2-2 p,m 

'02,2 

2oj + 2(U - n s ) 

I5/i s aVi s 2 2 

64a|(l-e|) 3 / 2 U+ ' 

02,1 

2 lo -|-i2 — f?g 

15/i S a 2 si s cis 2 ■/-, . 

16a3(l- e |)3/ 2 eSl ( 1 + ^ 

02,0 

2 uj 

15/isa 2 (3s 2 is — 2) 2 2 - 
64a|(l-e 2 )3/2 6 S * 

0-2,0 

—2oj 

15/i s a 2 (3s 2 i s - 2) 2 2 . 
64a|(l-e 2 )3/2 6 S * 

0-2,1 

—2c j f? — f?g 

15/i s a 2 si s ci s 2 

- ^ - y;- —— e SZ( 1 — Cl) 

16a|(l — eg) 3 / 2 

0-2,2 

— 2oj + 2(U - Us) 

I5ii s a 2 s 2 i s 2 2 

64a|(i_ e |)3/2 e < cl) 

00,1 

f? — f?g 

3ms« 2 s*sc*s 2 . . 

8a|(l-e 2 )3/2 (2 + 3e 

00,2 

2(J2 - U s ) 

3/i s a 2 s 2 i s 3e2)s 2 i 

32a3(l-e|)3/2 (2 + 3e jS 


the GRIM5-S1 model. Our semi-analytical propagator is configurable and the third-body perturbations 
from the Moon and the Sun have been developed up to degree 4 and 3, respectively. All of the equations 
of motion have been formulated through the equinoctial elements, related to the Keplerian elements by 


a, 

e y = esin(cj + 12), 

Cx = ecos(u; + 12), 
i y = sin(*/2) sin(l2), 
i x = sin(i/2) cos(12), 
£ = M T w T 12 , 


(46) 


which are suitable for all considered dynamical configurations. The variational system (i.e., the equations 
of motion for the state and tangent vector) are then propagated using a fixed step size (set to 1 day) 
Runge-Kutta 6-th order integration algorithm. 

To produce the different maps of the atlas, only the initial eccentricity and inclination were varried, 
distributed in a regular grid of 320 x 115 initial conditions for the four semi-major axes of Fig. [3] The 
computation of the FLIs on a grid of initial conditions allows a clear distinction, in short CPU time, of 


Lega et al. 

2010 

Guzzo 

CM 

o 

o 

CM 


Guzzo 2002). The number of initial conditions chosen 


was a good balance between CPU time and the final pictures offered by the resolution. Concerning the 
simulations, after a calibration procedure, we decided to present the results of the FLI analysis after 530 
years of propagation. By increasing the iteration time, the basic features of the FLI maps are not changed, 
but small higher order resonances can be detected in a proper resolution. Our chosen timescale may seem 
prohibitive and might also not be the best trade-off between clear distinction of the separation of nearby 
orbits (if ever) and CPU time. However, as showed by Fig. 15 a 250 years propagation may leads to partially 


erroneous conclusions concerning the stability in MEO. Moreover, for weakly hyperbolic orbits (usually for 
moderately eccentric orbits), the time necessary to detect the divergence is longer. Consequently, there is 
no risk about the conclusions by presenting the results after this propagation time, at the cost of somewhat 
more CPU time (T odorovic et al.(2 008). For the zoomed-in portion, the resolution and propagation duration 
have been increased to get finer details, structures, and more contrast in the maps. Unless explicitly stated, 
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Fig. 14: Resonance centers and widths for the special class of inclination-dependent-only lunisolar resonances 
(q.v. Hughes||1980), plotted both with and without the solar contribution. 


all others parameters (initial phases, initial epoch, initial tangent vector) have been fixed for each map. 
Concerning the initial tangent vector, we used the same for the whole map, a fixed normed vector orthogonal 
to the flow, /, in Eq. |36[ The robustness of our results with respect to the choice of wq was tested by 
computing the maximal Lyapunov exponent. The Lyapunov exponent, which is first of all a time average, 
is a property of the orbit, independent of the initial point of that orbit. By drawing the mLCE maps, we 
observe a nice agreement with the FLI maps in terms of the detected structures (chaotic or stable). The 
only (obvious) disagreement was the contrast in the maps, easily explained: Lyapunov exponents are slow 
to converge. This argument is strong enough to attest the relaxation that we can have with respect of wo- 
Finally, the final values of the FLIs are coded, and projected into the plane of i—e initial conditions, using 
a color palette ranging from black to yellow. Darker points correspond to stable orbits, while lighter colors 
correspond to chaotic orbits. The FLIs of orbits whose propagation stop before reaching the final time, 
tj, have not been considered in the maps: they appear as a white to avoid the introduction of spurious 
structures^] This is the reason why some maps seems to be perforated. Re-entry orbits were in practice 
declared if there exist a time t <tf such that o(l — e) < 120 km. 
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